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1.  Introduction 

Consider  the  statistical  model  A/--=fX,A,{P*:0€n  })  and  event  C€A  and  suppose  we 
wish  to  estimate  P,(C)  based  on  a;  sample  X~ (xu  .  .  .  ,x,)  from  M.  The  typical  approach  to 
this  problem  is  to  select  a  probability  measure  Qx  on  A,  dependent  on  the  observed  data,  and 
then  quote  Qx{  C )  as  the  estimate. 

For  example  the  nonparametric  estimator  of  Pi(C)  is 

.  (i) 

ni-u 

where  Ic  is  the  indicator  function  of  C;  i.e.  the  empirical  probability  content  of  C.  For  a 
sufficiently  broad  class  {Pt:8€Cl  }  this  estimator  is  known  to  be  UMVU. 

If  m(X)  is  a  complete  minimal  sufficient  statistic  for  {.P|;0€fi  }  then 

<tf(C)-£(<tf(C):m(JO|  (2) 

is  UMVU.  Clearly  is  a  probability  measure  on  A  as  it  is  formed  by  mixing  probability  dis¬ 
tributions. 

Perhaps  the  most  commonly  used  method  of  obtaining  an  estimate  is  to  choose  some 
estimator  8(X)  of  8;  e.g.  the  MLE,  and  then  quote 

<?/(<?)-*WC)  (3) 

as  the  estimate. 

In  a  Bayesian  context,  or  perhaps  just  as  a  method  of  generating  a  plausible  estimate, 

a  prior  for  9  gives  rise  to  a  posterior  Px  for  6  which  in  turn  induces  a  distribution  for  Pt(  C). 

Then  the  minimum  Bayes  risk  estimator  with  respect  to  squared  error  loss  is  given  by 

Qg(C)-f  P,(C)dPx(9)  (4) 

o 

the  expected  value  of  P ,  with  respect  to  the  posterior.  Again  Q§  is  a  probability  measure  on  A 
as  it  is  a  mixture  of  probability  measures. 

Other  strategics  could  also  be  devised  for  abtaining  estimates  but  we  will  restrict  our 
discussion  to  those  presented  above.  In.  all  of  these  approaches  we  note  that  the  choice  of  Qx 
docs  not  depend  on  C.  As  such  it  teems  more  appropriate  to  say  we  are  estimating  Pt  rather 


than  P,(C). 


One  way  of  inducing  al  least  some  dependence  on  C  is  via  the  joint  invariant  group  of 
the  model  M  and  the  event  C;  namely  the  class  G  of  those  1-1,  bimeasurable  g  :X  —*  X  satis¬ 
fying  gC*=  C  and  both  of  Ptg  and  P»g~x  are  in  H  }  for  all  0€fl .  If  X  also  satisfies 

topological  requirements  then  it  makes  sense  to  require  that  g  also  preserve  this  structure;  e  g. 
if  X—Rr  then  we  require  g  to  be  a  diffeomorphism  (1-1,  onto  and  infinitely  differentiable 
both  ways).  If  x~Pj  then  gx~P#0*1  and  P#(  C)*=P|(  g~  xC)**P,g~ '( C).  Thus  the  estimate 
of  Pf(C)  should  satisfy  Qx{C)’^‘QtX(C)i  i.e.  our  estimate  should  be  the  same  whether  we 
observe  X  or  gX.  As  we  shall  see,  this  criterion  leads  to  some  restriction  in  the  class  of  possi¬ 
ble  estimators  for  the  problem  we  consider  in  the  succeeding  sections. 

2.  Circular  Error  Probabilities  for  the  Bivariate  Normal 

Suppose  that  x~N2(/i,E)  and  C4«={x:  x'x<  k2}.  Thus  x  could  give  the  coordinates 
of  the  hitting  point  for  some  projectile  aimed  at  buliseye  0  and  we  wish  to  estimate  the  proba¬ 
bility  of  coming  within  k  of  o  as  an  assessment  of  the  accuracy  of  the  targeting  procedure. 

Even  when  (;i,E)  is  known  the  problem  of  calculating  P(*x)(CV)  ia  significant.  For 
various  tabulations  and  results  related  to  this  problem  see,  for  example,  Grad  and  Solo- 
mo^  1955),  HarterflOdO),  Lowe(1980),  Groenewoud  et  al.(1967)  and  Govindarajulu(1983). 
For  further  problems  involving  probability  calculations  related  to  targeting  problems  see,  for 
example,  Solomon(  1953)  and  Guenther  and  Terragno(l904). 

The  problem  we  are  concerned  with  here  is  to  estimate  P(„,e)((7*)  based  on  a  sample 
A'w^x,,  .  .  .  ,x„)  from  the  N2(/i,E)  distribution  where  and  Eg/?2*2  positive  definite 

are  unknown.  The  relevant  invariant  group  for  this  problem  is  0(2),  the  group  of  orthogonal 
transformations  on  R2.  In  particular  we  discuss  algorithms  for  the  evaluation  of  estimators  of 
the  form  (2),  (3)  and  (4)  which  satisfy  the  invariance  requirement.  Further  we  present  some 
Monte  Carlo  results  which  give  some  indication  as  to  the  bthaviour  of  the  estimates  in 
repeated  sampling  and  as  such  present  additional  information  for  the  investigator  who  might 
be  faced  with  choosing  amongst  them. 
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The  estimators  appropriate  to  the  situations  when  (/j,E)  is  restricted;  e.g.  requiring 
that  E“<r2/  and  a  >0  unknown,  can  typically  be  obtained  by  making  obvious  adjustments  to 
our  algorithms  for  the  most  general  case.  The  computer  programs  for  the  evaluation  of  the 
estimators  and  the  simulation  were  written  in  Fortran  77  and  are  available  from  the  authors. 

All  the  calculations  discussed  in  this  paper  were  carried  out  on  the  FDP  11-70  in  the  Depart¬ 
ment  of  Statistics,  University  of  Toronto. 

The  dual  of  the  problem  addressed  here  is  to  specify  PommP(^,r)[Cls)  and  then  based 
on  the  data  X  estimate  k.  When  po**».5  the  value  of  k  is  referred  to  as  the  circular  probable 
error.  This  problem  is  discussed  in  Blischke  and  Ralpin(  1966). 

8.  The  Standard  Estimator 

By  the  standard  estimator  we  mean 

-.)-«<*>  <5' 

where  xs=n*‘J]x,  -ad  5x-=-(X-xl')(X-xl')';  i.e.  we  have  replaced  n  and  £  by  their  UMVU 
•- i 

estimators.  This  estimate  is  clearly  invariant  under  0(2).  The  tabulations  mentioned  earlier  are 
available  for  the  calculation  of  (5).  This  approach,  however,  requires  interpolation  and  is  not 
appropriate  for  extensive  Monte  Carlo  work.  We  discuss  two  approaches  to  the  computer 
evaluation  of  this  estimator. 

First  we  write  S’^=»SASA'  where  is  the  unique  lower  triangular  matrix  with 

positive  diagonal  elements  satisfying  this  equation.  Then  if  s~Af2(o,/)  we  can  write  (5)  as 

Hx+(n-l)->*S*s||  <  k)  =P (o.;)(  £«?(*+»,)*  <  ("-1  )k‘ )  (6) 

•  —1 

where  b=(  n-  l)'^^.^1  x  and  QDQ  is  the  spectral  decomposition  of  S^S^. 

An  algorithm  for  calculating  (6)  is  obtained  by  using  the  fact  that  i i(  *3~JV(0,i),  sta¬ 


tistically  independent  and  thus  we  can  write  (6)  as 


•»'*  l  V  *-%  i  "V  -»r±  »y_  4  — <_  4VJ.^: 


j 
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*> 

►  ' 
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-*a+^»(.-i)>^* 

/  [*(«(-))-*(/(1))^(,)rfs  (7) 

-*r<Jl(-  0,/!>* 

where  <t  is  the  distribution  function  for  the  N(0,1),  is  the  density, 
and  u(t),  l(z)  equal 

-  i|d:  df1  [(n-  l)ia-  d|  (*»+&a)*l l/*  (8) 

respectively.  An  efficient  algorithm  is  then  obtained  by  using  a  packaged  routine  for  the 
evaluation  of  4>,  e.g.  IMSL,  and  performing  the  integration  using  a  Gauss-Legend  re  rule. 

A  disadvantage  of  the  above  approach  arises  when  we  are  interested  in  the  higher 
dimensional  analogs  of  this  problem  as  the  computation  becomes  progressively  more  compli¬ 
cated.  A  more  efficient  approach,  and  it  is  the  one  we  have  adopted,  is  based  on  an  adaptation 
of  an  algorithm  due  to  Sheil  and  0'Muircheartaigb(  1977)  which  is  in  turn  based  on  results  due 
to  Ruben(1962)  concerning  the  evaluation  of  the  distribution  of  ( *+b) 'D  ( s+b)  where 
s~A/^(o,/),  bG-R*  and  D&RpX*  is  diagonal  with  nonnegative  diagonal  elements.  The  result 
gives  a  series  representation  for  the  distribution  function  of  this  quantity  and  thus  aUo  a  aeries 
representation  for  (6).  We  controlled  the  accuracy  in  our  calculation  by  stopping  the  summa¬ 
tion  when  the  contribution  of  the  remaining  terms  was  less  than  I0~n.  For  p=»  2,  n»  20  an 
evaluation  of  the  estimate  takes  approximately  .1  seconds  of  CPU  time. 

<L  The  UMVU  Estimator 

We  have  that  (x,  Sy)  is  a  complete  minimal  sufficient  statistic  and  thus  the  UMVU 
estimator  is  given  by 

«(T(C*)-£|^(C*)  :  x,  5x]-£(/ct(x,)  :  x,  Sx\.  (9) 

This  expectation  can  be  evaluated  using  a  result  due  to  Laurent(i957);  namely  the  conditional 
distribution  of  X)  given  (x,  Sy)  has  density  proportional  to 

[l-(l-l/»)-,(xl-x),5jf,(x,-S)]‘-‘Ws  (1C) 

where  x(  is  constrained  so  that  the  term  in  the  brackets  is  positive.  If  we  transform  X—*QX 
where  QGO(2)  we  see  that  (10)  is  unchanged  and  thus  this  estimator  is  invariant  under  0(2). 
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Laurent  proposed  chat  a  calculation  such  as  (0)  could  be  carried  out  by  a  double 
numerical  integration.  In  fact  we  can  simplify  the  calculation  of  (9)  subs  tan  dally.  First  we 
moke  the  transformation  x,— »t  where  t=  Vl-  1/n  The  density  of  t  is  uien  propor¬ 

tional  to 

(11) 

Then  in  a  particular  quadrant  of  R3  we  make  the  transformation  t-*v  where 

ijn(t,).  This  transformation  has  Jacobian  ( l/2)2vf  and  thus  the  density  of  v,  con¬ 

ditional  on  a  quadrant,  is  proportional  to 

(12) 

From  this  we  conclude  that  (vi,o2,l-  »i-  v2)  is  distributed  D irictilet(  1  /2, l/2,( n—  3) /2) . 

The  estimate  can  then  be  expressed  os 

T  £  />{(*- WI^  «,»<*, 

i 

(Ij+n/i-  l/n  sjidtVj^+v/l-  l/n  Saadjvj^1)2^  fcs}  (13) 

where  5A"=(«,y)  and  P  refers  to  the  distribution  of  (®j,v-).  We  now  make  the  transformation 
( ®i i ®s) -*(“,«)  where  u»vi/(l-v2)  and  »«=»«2,  This  transformation  has  Jacobian  i-v  and  thus 
the  joint  density  of  (u.v)  is  proportional  to 

«',/a(  1-  */»(  1-  ( 14) 

Therefore  *~Beta(  1/2, ( n-  3)  /2)  is  statistically  independent  of  »~Beta(l/2,(n-2)/2).  Then 
denoting  the  Beta(p,q)  distribution  and  density  functions  by  B[  '  |  p,q)  and  6(  ■  |  p,q) 
respectively  (13)  can  be  written  as 
, 

T  £  /iB(tio(v)|l/2,(n-3)/2)-B(«I(v)|l/2,(n-3)/2)|4(v|l/2,(n-2)/2)</v  (15) 

■*  Ma”±  >  *i 

where  putting  r2,  r1  equal  to 

(d2(*  (Jjl-  *2«n)±  («n  +«2l  )l^)/v/l- l/n  *11*22  (19) 

respectively,  we  have  that 

»(  =  |mm(ma*(r.-,0),i)]2  (17) 


and  putting  <2,  s,  equal  to 


(-  4(w)±(6i(«)-4oe(v))l^]/2«Vrl-  v 

respectively  where 

‘"(1-  l/n)(«?i  +4ji) 
i(e)—2\/l-  l/n  rfi(*iSn+*a4ai+'A~  1/n 

*( »)»( *a+v^l-  l/w  4aa^awl/1)2"K*?-  O 

then  we  have  that 

«<(  • )  —  t  mm  ( ma  at  ( s,  ,0) ,  1 )  ] 2. 

To  evaluate  (15)  we  use  an  IMSL  subroutine  for  B{  •  |l/2,(n-3)/2)  and  then  use 
Gauss-Legendre  integration.  The  efficiency  of  the  integration  is  improved  by  using  an  IMSL 
routine  for  the  inverse  of  the  Beta(  l/2,(  n- 2)/2)  distribution  function  to  find  the  point  vsup 
such  that  B(vsup|l/2,(n- 2) /2)-«  .999999  and  then  using  min(oa  ,vsup)  as  the  upper  limit  in 
the  integration.  This  ensures  that  the  Gauss  points  are  concentrated  where  the  probability  lies. 
We  further  improved  the  efficiency  of  the  integration  by  making  the  transformation  v—w 
where  ui  —  v1^  so  that  w  has  density  proportional  to  (l- This  transformation 
removes  the  singularity  atO  which  b(  •  |l/2,(  n-  2)/2)  possesses. 

The  accuracy  of  the  calculation  is  controlled  by  dividing  the  interval  of  integration 
into  subintervals  of  equal  length  and  using  a  Gauss  rule  of  the  same  order  in  each.  The  pro* 
gram  allows  the  order  of  the  rule  to  vary  from  1  to  20  and  for  as  many  subdivisons  as  desired. 
Thus  arbitrary  accuracy  can  be  achieved  with  the  tradeoff  being  computing  time.  For  n=»  20, 
using  4  subdivisions  and  a  Gauss(lO)  rule,  stability  was  achieved  in  the  fifth  decimal  place  and 
took  about  1.5  seconds  of  CPU  time. 

6.  Bayes  Estimate 

There  are  of  course  many  different  Bayes  estimates  as  there,  are  many  different 
choices  for  the  prior  distribution  of  0  and  for  the  loss  function.  Here  we  will  use  mean-square 
error  and  choose  the  prior  distribution  for  (n, E)  to  be  Jeffrey’s  prior;  sec  for  example  Box  and 
Tiao(l973). 

An  alternative  approach  to  this  problem,  which  also  leads  to  the  Bayes  estimate 


the  itf'Ctural  model  for  the 


associated  with  Jeffrey's  prior,  is  to  use 
multivariate  normal  model  based  on  the  affine  group;  i.e.  the  group 

G-{la,C]  :  a€i?3,  Cert2*3,  *<(  C)^  0},  (»i,C,,|(«a,Ca]  —  and  we  represent 

x~Ndi*£)  as  x-[/i,r]s-p+rs  where  T  satisfies  E-fT’  and  s-N^o,/).  Utis  approach  is 
discussed  in  Fraser(l979)  where  the  positive  affine  group  is  used;  i.e.  we  also  require 
det(C)  >0.  Use  of  the  full  affine  group  requires  only  minor  adjustment  to  the  analysis 
presented  in  Frascr(1979)  and  it  provides  a  convenient  framework  for  obtaining  our  results. 

The  structural  model  leads  to  the  following  relations 

r— cxc-‘ 

where  E-IT',  TeR3*  with  det(T)^  0,  SX-CXCX'  with  CxeR3y2  calculated  as  described  in 
Fraser^  1979) ,  ,n~lI)  statistically  independent  of  C  which  has  density  as  described  in 

Fraser<l979)  with  the  adjustment  that  the  density  is  multiplied  by  2-3  and  is  now  a  function 
on  [C  :CeR3*,  det(C)^O). 

As  is  well-known,  the  Bayes  estimate  with  respect  to  mean-square  error,  is  given  by 
the  mean-value,  assuming  it  exists,  of  the  marginal  posterior  distribution  of  the  quantity  to  be 
estimated.  Thus  we  wish  to  calculate 

Denoting  the  joint  posterior  distribution  of  (pj)  by  Px  and  using  E^IT'  we  can  write  (22) 

as 

/*WllM+r»||< 

This  expectation  can  be  evaluated  by  using  the  relations  (21)  and  the  joint  distribution  P  for 
(s,  C)  to  obtain 

j  ./)(  l^v^T  CyMI—  k)dP  (s,(7) 

where  t-Vl+l/n  C-'(s-I).  From  Frascr(1979)  we  have  thatthas  density  function  given  by 
where  d,=2ir/^/T(//2);  i.e.  t  has  a  canonical  bivariate  Student  (n-1)  distribution.  Using  the 


Gram-Schmidt  decomposition  on  the  rows  of  C*  we  obtain  Cx~S&Q  where  Sa  is  as  before 
and  Q€0{2).  We  note  that  the  distribution  of  t  is  invariant  under  orthogonal  transformations. 
Then  by  the  theorem  of  total  probability  (24)  can  be  written  as 

P((F1+\/l+l/«»  *n<i)!+(*j+s/m7n"  aia*i+'/l+l/n  ^aa^)3  —  ^S) 
where  P  refers  to  the  distribution  of  t 

We  now  make  the  transformation  ( f r t *a)  — 1 ’(u»v)  where  «*'/n-l  ft/V^-b^a 
v~Vn^~2  tj.  The  density  of  (u,v)  is  proportional  to 

and  thus  u~Student(n-l)  statistically  independent  of  v~Student(n-2).  Therefore,  denoting 
the  Student X)  distribution  and  density  functions  by  G(|X),  f(-|X)  respectively,  we  have  that 
(27)  can  be  written  as 

•« 

/  [C(«2(«)  |«-l)-<7(«i(v)  |n-l)]p(u  |n-  2)  iv 

«i 

where  v2,tq  equal 

y/n-2  [  *a4ud:  *  V *?i  +*ai  ]  /n/1+1/"  »n4aa 

respectively,  equal 

[l+„2/(n-2)]-l^l-  h1(V)=fcCh?(  v)-4«,c1(v))1^]/2<al 

respectively  and 

•  i  =  (l+l/n)(a?i+«|i ) 

b  1(v)=2V’H-l/n  [(F|Su+*a2vT+l/n  [(*i*u+F2»2i)  Wl  +  l/n  2  ] 

e,(v)=(  Fj+v/l  +  l/n  Sojv/v'n-2  )s-rz?-k' 

The  calculation  is  then  carried  out  using  an  DvlSL  subroutine  for  the  Student  distribu* 
tion  function  and  Gauss-Legcndre  rules  for  the  integration.  The  accuracy  was  controlled  by 
subdividing  the  interval  (vl,v2)  into  subintcrvals  of  equal  length,  carrying  out  the  numerical 
integration  within  cr.ch  subinterval  using  a  Gauss  rule  of  the  same  order  and  controlling  the 
number  of  points  in  the  rule.  If  V|<0<v2  then  each  of  the  intervals  (v,,0)  and  (0,t2)  was 
subdivided  into  the  same  number  of  subintcrvals  of  equal  length.  This  is  to  ensure  that  the 
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mode  0  of  the  Student  density  serves  as  an  endpoint  for  the  integration  as  this  improves  the 
efficiency  of  the  calculation.  A  further  improvement  was  made  by  requiring  that  jr,-|  be  no 
greater  than  <.oooooi( n_  2)  and  this  point  was  obtained  from  an  IMSL  subroutine  for  the  inverse 
of  a  Student  distribution  function.  For  n»*  20,  using  4  subdivisions  and  a  Gauss(lO)  rule,  star 
bility  was  obtained  in  the  fifth  decimal  place.  This  calculation  took  approximately  .8  seconds  of 
CPU  time. 

As  discussed  above,  the  estimator  (28)  is  ootained  using  Jeffrey's  prior  and  the  result 
is  also  obtained  from  the  structural  model  using  the  affine  group.  In  effect  Jeffrey’s  prior 
results  as  the  marginalization  of  the  right  Haar  prior  on  this  grou^.  As  is  well-known  other 
groups  can  be  used  to  parametrize  the  multivariate  normal  model  and  their  right  Haar  priors 
give  rise  to  different  priors  for  the  full  parameter  (#i,£).  For  example,  the  affine  lower  triangu¬ 
lar  group  leads  to  the  estimate 

P(||*+v'I+T7"‘  S*t\\<  k)  (32) 

where  P  refers  to  the  distribution  of  t  which  has  density  proportional  to 

(33) 

by  results  in  Fraser(  1979) .  Thus  we  see  we  will  obtain  a  different  estimate  of  a  form  similar  to 
(28)  in  this  case. 

The  invariance  considerations  lead,  however,  to  the  choice  of  the  estimator  (28).  For 
if  we  transform  X  to  QX  where  Q€0(2)  then  (x.Cy)  transforms  to  (Q*,Q(?x)  and  from  (24) 
we  see  that  the  estimator  is  invariant  under  0(2).  This  invariance  property  does  not  hold, 
however,  for  the  estimator  based  on  the  affine  lower  triangular  group.  For  example,  if  Sx—  I 
then  5^=1,  QSxQf  —  I  and  (33)  is  not  invariant  with  respect  to  0(2)  which  proves  the  non- 
invariancc  of  the  estimate  in  this  case. 

#.  The  Monte  Carlo  Study 

To  study  analytically  the  repeated  sampling  behaviour  of  the  estimates  we  have  dis¬ 
cussed  presents  a  difficult  problem.  Accordingly  a  simulation  study  was  carried  out  to  see  how 
effective  the  estimators  are  and  tc  assess  their  relative  merits. 
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The  performance  of  the  estimators  was  considered  for  four  sets  of  parameter  values 

(0 

(■)  —  S— / 

( m)  E=.511'+.5/ 

(»)  p  — 1,  E— .5ll'+.57. 

For  each  parameter  set  we  calculated  k  such  that  P(^,q(  C* )  — .5  using  the  algorithm  of  Shcil 
and  0’Muircheartaigh(l977).  Our  estimates  were  then  always  of  the  true  value  .5.  For  each 
parameter  set  we  considered  the  estimation  problem  for  sample  sixes  n~  10  and  n1—  20. 

For  a  given  parameter  (p,E)  and  sample  sixe  n  we  do  not  need  to  generate  the  full 
sample  Xa>(x1>  .  .  .  ,x.)  from  N^^,Z)  to  calculate  the  estimates.  For  we  need  only  generate 
(ff.Sx)  where  ST — IV^/j^'E)  statistically  independent  of  S* — Wj(E,n-l).  To  do  this  we  used 
the  following  relations 

7— ii+n-  l^EAa 

Sjf“EASA5A'EA'  (34) 

where  a— N^O,/)  statistically  independent  of  Sa’“,(4»/)  'where  *?i  ‘~Chi-square(n-l), 

#22 ~Chi-square( n-2) ,  #31- — N(0,1) ,  #,a*«0  and  «u,  #2S>  *ai  statistically  independent.  The 
N{0,1)  variablea  were  generated  using  the  Box-Muller  method;  namely 
#“*(-2foj(«1))l^cos(2ir«a)  where  «i,  we  statistically  independent  and  diatributed  U(0,1). 

The  chi-square  variables  were  generated  using  a  method  due  to  Cheng  and  Fcast(  1979)  and  in 
fact  we  used  the  program  RGKM3  as  it  is  listed  in  Bratley,  Fox  and  Schrage(1983) . 

The  uniform  random  variates  needed  for  the  generation  of  the  normals  and  chi- 
squares  were  obtained  using  the  routine  due  to  Schrage(  1979) .  To  decrease  possible  effects 
due  to  noUquite-randomness  we  first  filled  a  table  with  uniform  values  .  Then  each  time  a 
value  was  required  we  generated  a  random  address  ,thc  contents  of  which  becomes  the  gen¬ 
erated  value,  and  replaced  it  in  the  table  by  a  newly  generated  value. 

For  any  given  (x,5^),  generated  as  above,  we  computed  oil  three  estimators;  i.e  we 
used  common  random  numbers.  Accordingly  the  estimators  share  equally  in  any  effects  due  to 


-3  :j  ** ->  rj>  ?  i- 
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dcficicncics  in  the  generators.  Further,  as  we  will  aeer  this  technique  substantially  improved 
the  efficiency  of  the  Monte  Carlo  study. 

For  a  given  sample  site  and  parameter  set  we  generated  nrep  values  of  (£,5*)  where 
nrep  depended  on  the  sample  site  n.  Then  for  each  estimator,  denoted  genetically  for  con¬ 
venience  by  x,  we  estimated 

Ae(*)  —  E(t) 

A/S£(*) -£(*-.&) 2 

by 

-5)3 

nrtP  <- 1 

respectively.  The  standard  errors  of  these  estimates  are  given  by 

SD  ( Av(  x ) )  =  (£(  x2)  -  ( £(  s ) ) *] ‘'V^F 
SD(MSE{  z  ))  —  [£(*-  .5)*-  (2?(r-  .S)2)2]  '^/Vrirep 
respectively  and  they  in  turn  are  estimated  by 

&>(viv(z))  — [  Av(*))2/nrep(nrep-l)l1'* 

t— I 

£d{MSE{x)  —  [  Xf((x>~  &)J-  MSE(x))i/nrtp(nnp-  1*)1 
f-i 


(35) 


(36) 


(37) 


(38) 


respectively.  We  then  tested  the  null  hypothesis  that  x  is  unbiased  for  .5  via  a  t-test  using  the 


statistic 

**-(  At>(i)- .5)/§D(Av{x))  (39) 

i.e  we  compare  this  value  with  the  N(0,l)  distribution  by  computing  the  observed  level  of 
significance  P(|Z  |  >|z  |)  where  Z—  N(0,1). 

Ibe  primary  purpose  of  the  study  was  to  compare  MSE(x)  with  MSE(y)  for  estima¬ 
tors  x  and  y.  When  x  is  the  nonparamclric  estimator  we  have  that  MSE(x)>=»  ,25/nrep  and  the 


test  statistic  takes  the  form 
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*-( lClSE(y)- MSE{  z))/SD(  tiSE(  y ))  (40) 

When  x  and  y  are  two  of  the  estimators  we  discussed  in  the  preceding  sections  the  test  statistic 
takes  the  form 

i - ( A lSE(y)-  fiSE{ x))/SD(  fiSE (y)-tiSE(z))  (41) 

where 

^(Adr5E(|r)-A/5£(I))-(^D^(A/5£(y))+5f)2(i^S£(J))-2<^,o»(A3'5£(>),MS^(*))],^  (42) 

is  the  estimate  of 

SD(MSE(y)~  K(SE(x))**[(  Var(y-  .5)S+Var(x-  .5)1- 2Cot>((y-  .5)*,(x-  .5)2))/nrep] (43) 
sad  where 

(5’o»(A5r5£(y),AV5£(x))»"^((y/-  .S)^*^-  .h)3-  tkSE{y)MSE(x))/nrcp{nrcp-  1).  (44) 


The  covariance  term  is  required  in  (43)  because  we  have  used  common  random  numbers. 
From  this  we  see  where  die  gain  in  efficiency  was  obtained  as  in  all  cases  MSE(x)  and 
tiSE(Y)  were  positively  correlated  and  this  reduced  SD(trfSE[y)-  MSE[x))  substantially. 
The  estimated  correlation  between  these  two  quantities  ranged  from  .460  to  .096  and  in  most 
cases  was  greater  than  .700. 

The  number  of  replications  for  each  sample  site  and  parameter  set  was  determined  by 
first  performing  a  trial  run  of  100  for  all  cases  and  calculating  the  estimates  we  have  just 
described.  The  primary  determinant  of  sampling  variability  turned  out  to  be  the  sample  size. 
We  then  estimated  an  upper  bound  for  $D(MSE(x))  for  all  estimators  over  all  parameter  sets 
within  a  sample  size.  On  the  basis  of  this  information  we  chose  nrep  so  that  when  n*=  10  the 
half-length  of  a  ,95-confidrnce  interval  fr>r  MSE(x)  would  be  less  than  .001  and  when  n*=  20 
so  that  a  ,95-confidcnce  interval  would  have  half-length  less  than  ,00025.  With  some  extra 
margin  for  safety  this  lead  to  choosing  nrep—  15000  when  n—  10  and  nrep—  10000  when 
n*«  20.  The  results  we  obtained  tend  to  confirm  our  expectations.  These  choices  gave  con¬ 
clusive  results  for  the  comparisons  amongst  the  MSE(x)  because  of  the  use  of  common  ran¬ 
dom  numbers.  Further  these  values  of  nrep  gave  that  the  half-length  of  a  .005-confidencc 
interval  for  Av(x)  is  less  than  .005. 
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For  each  estimator  the  accuracy  of  the  calculation  was  controlled  so  that  the  error  was 
less  than  5X10'*  in  each  evaluation  of  the  estimate;  i.c.  if  £  denotes  the  computed  value  of 
the  estimate  and  z  *  denotes  the  actual  value  of  the  estimate  then 

5X10-*  (45) 

Thus  the  absolute  error  in  Av(z)  is  less  than  5X10"*  and  since 

X,*-.S)*-(i~.5)3|<  KOa-i)2)+2|**-i  |  (46) 

we  have  that  the  absolute  error  in  KtSE{z)  is  less  than  2)40**.  All  calculations  were  done  in 
double  precision. 
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PinmtMn  |  Estimator 

Av(x)  | 

(i) 

Laurent 

0.501630 

Standard 

0.507041 

Bt,tt» 

0.422713 

(«) 

Laurent 

0.409720 

Standard 

0.465583 

Bajtt 

0.443320 

(ill) 

Lanrtnt 

0.501258 

Standard 

0.505730 

Bajrtt 

0.422180 

(W) 

Laurent 

0.400300 

Standard 

0.401079 

Btjrn 


0.434217 


Tablt  1 


Sam  pit  tlit 

-  10 

$D{Av{x)) 

<*A 

tASE{  x ) 

SD  (  A/S£(  X ) ) 

au 

0  00004870 

0.085 

0.013505 

0.00016076 

0.0 

0.00088234 

0.0 

0.011740 

0.00014092 

0.0 

0.00077622 

0.0 

0.015011 

0.00013103 

0.0 

0.00102002 

0.787 

0.016806 

0.00018432 

0.00101844 

0.001 

0.016669 

0.00018072 

0.00087711 

0.0 

0.014752 

0.00016769 

0.00002120 

0.174 

0.127320 

0.00016824 

0.00086404 

0.0 

0.011231 

0.00013010 

0.00076096 

0.0 

0.014879 

0.00012916 

0.00100313 

0.640 

0.015093 

0.00018800 

0.00100830 

0.0 

0.016329 

0.00018626 

0.00058902 

0.0 

0.016182 

0.00017134 

Tablt  2 

Samplt  tii««  10 

ParamtUn 

Comparison  | 

U 

$D(MSE{y)- MSE(x}) 

(») 

Laurtnt  t*  Standard 

74.8 

0.00002383 

Laurtnt  ti  Bajrtt 

-10.2 

0.00014714 

Standard  vi  Bayti 

•23.2 

0.00014008 

(a) 

Laurtnt  ti  Standard 

23.2 

0.00001572 

Laura  ntta  Bayta 

73.5 

0.00011812 

Standard  rl  Bajrtt 

77.5 

0.07010625 

(Ul) 

Lauraat  «i  Standard 

88.5 

0.00002257 

Laurtnt  ti  Bajrta 

•14.7 

0.00014610 

Standard  ra  Bajrta 

•286 

0.00013780 

(M 

Laurtnt  ti  Standard 

11.4 

0.00002084 

Laurtnt  vi  Bayta 

•  8.33 

0.00013084 

Standard  ts  B nytt 


-7.47 


0.00011428 
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|  TO.lt  3 

Simple  slit*  20 

Parameters 

Estimator 

Av(x) 

§D[Av(z)) 

u 

MSE[z) 

$D(MSE(x)) 

(l) 

Laurent 

0.S00707 

0.00070051 

0.882 

0.006302 

0.00009100 

0.0 

Standard 

0.50  4093 

0.00077234 

0.0 

0.005989 

0.00008087 

0.0 

Bayes 

0.400413 

0.00072084 

0.0 

0.008783 

0.00008111 

0.0 

(») 

Lanrent 

0.500332 

C. 00087 090 

0.529 

0.007584 

0.00010944 

0.0 

Standard 

0.408784 

0.00087101 

0.181 

0.007587 

0.00010904 

0J> 

Bayea 

0.473130 

0.00080048 

0.0 

0.007225 

0.00009973 

0.0 

(IU) 

Laurent 

0.500007 

0.00077780 

0.240 

0.008050 

0.00009016 

0.0 

Standard 

0.503080 

0.00075381 

0.0 

0.005807 

0.00008533 

0.0 

Bayes 

0.400345 

0.00070728 

0.0 

0.000575 

0.00007921 

0.0 

(»*) 

Laurent 

0.500133 

0.00085500 

0.89 

0.007309 

0.00010730 

0.0 

Standard 

0.405074 

0.00085777 

0.0 

0.007376 

0.00010738 

0.0 

Bayes 

0.485767 

0.00081321 

0.0 

0.007784 

0.00010672 

0.0 

Tibia  4 


Sun  pit  ill*—  30 


Parameters  j 

|  CompAritoD  | 

■ 

$D(MSE{y)-MSE{z)) 

(1) 

Laurent  vs  Standard 

52.8 

0.00000765 

Laurent  vi  Bayes 

-6.73 

0.00006470 

Staadwd  ti  Bayes 

-11.65 

0.00006607 

(U) 

Laurent  vi  Standard 

-0.74 

0.00000448 

Laariat  vi  Bayes 

7.44 

0.00004820 

Standard  vt  Bayea 

8.06 

0.00004400 

(111) 

Laurent  vt  Standard 

63.48 

0.00000859 

Laurent  vt  Bnjrta 

-8.22 

0.00006380 

Standard  vi  Bayes 

-13.56 

0.00006470 

(i») 

Laurent  vi  Standard 

-8.32 

0.00000800 

— 

Laurent  vt  Biyes 

-7.53 

0.00005991 

Standard  t o  Bayes 


-7.70 


0.00005308 


In  every  ease  it  turns  out  that  wc  have  no  evidence  against  the  hypothesis  that  the 
Laurent  estimator  is  unbiased  for  .5  and  this  is  os  theory  predicts.  We  see  that  in  every  case 
except  n**  20, (ii)  we  reject  the  hypothesis  that  the  Standard  estimator  is  unbiased  for  .5.  We 
note,  however,  that  the  bias  in  this  estimator  is  quite  small  in  every  cose  with  the  largest  esti¬ 
mate  of  the  bias  being  about  .000  and  the  bias  decreases  as  n  increases.  In  every  case  we  reject 
the  hypothesis  that  the  Bayes  estimate  is  unbiased  with  the  smallest  estimate  of  its  bias  being 
about  .027.  The  bias  decreases  as  n  increases  and  can  be  severe  for  small  sample  sites. 

In  every  case  we  reject  the  hypothesis  that  the  mean-square  error  of  the  estimator 
included  in  the  study  was  equal  <o  that  of  the  nonparametric  estimator.  The  Laurent,  Stan¬ 
dard,  and  Bayes  estimators  would  all  appear  to  be  substantial  improvements  over  the  non¬ 
parametric  estimator.  The  smallest  relative  efficiency  ,  as  measured  by  the  ratio  of  the  mean- 
square  errors,  of  an  estimator  to  the  nonparametric  estimator  was  154% 

We  now  compare  the  mean-square  errors  of  the  estimators  included  in  the  study.  We 
note  that  in  every  case  except  for  n— «  20,  (ii)  Laurent  versus  Standard,  we  categorically  reject 
the  hypothesis  that  the  mean-square  errors  are  equal.  For  (i),  in  both  sample  sizes,  we  have 
that  the  Standard  estimator  is  superior  to  Laurent's  which  is  in  turn  superior  to  the  Bayes  esti¬ 
mator.  For  (ii),  the  Bayes  estimator  is  superior  to  the  other  two  while  the  Standard  is  superior 
to  Laurent's,  when  n*»  10  and  they  are  equivalent  when  n»  20.  For  (iii),  we  have  the  same 
ranking  as  in  (i).  For  (iv).  Laurent’s  estimator  was  best  followed  by  the  Standard  which  in 
turn  was  better  than  the  Bayes  estimator  and  this  applied  for  both  sample  sizes. 

We  see  from  the  above  discussion  that  no  estimator  can  be  categorically  accepted  or 
rejected  as  the  best  or  worst  in  the  circumstances  we  considered.  On  the  other  hand,  when  tak¬ 
ing  account  of  both  bias  and  mean-square  error,  it  would  seem  that  the  Standard  estimator 
would  be  the  most  practical  choice.  In  fact  the  lowest  relative  efficiency  of  the  Standard  esti¬ 
mator  to  the  best  estimator,  when  it  was  not  best,  was  about  94%  The  lowest  relative 
efficiency  of  Laurent’s  estimator  to  the  best  was  about  87%  and  the  corresponding  value  for 
the  Bayes  estimator  was  about  75%  Perhaps  most  surprising  in  our  results  was  the  good  per- 
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formonce  of  the  Standard  estimator  relative  to  Laurent’s  estimator  given  that  the  latter 
possesses  on  optimality  property.  A  further  point  in  favour  of  the  Standard  estimator  is  given 
by  the  fact  that  a  much  more  cflicient  algorithm  is  available  for  its  evaluation  than  for  the 
other  two. 

7.  Conclusions 

This  paper  has  been  concerned  with  the  problem  of  estimating  circular  error  probabili¬ 
ties  when  we  require  that  the  estimator  be  invariant  under  the  invariant  group  of  the  circle. 
Three  competing  estimators  were  proposed  and  we  developed  efficient  methods  for  their 
evaluation.  A  Monte  Carlo  study  was  carried  out  to  provide  more  information  concerning  the 
relative  merits  of  the  estimators.  On  the  basis  of  this  study  and  the  relative  efficiencies  of  their 
algorithms  a  recommendation  can  be  made  that  the  Standard  estimator  ia  perhaps  the  most 
practically  useful  for  this  problem.  In  all  cases  the  estimators  were  substantially  better  than  the 
nonpar araetric  estimator  w  hen  we  are  assuming  bivariate  normality. 
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The  general  problem  we  are  concerned  with  here  is  the  estimation  of  P»(C)  whX^e 
C  is  some  fixed  event  and  Pg  is  unknown  in  some  class.  The  Various  available  es-\ 
timation  procedures  seem  to  involve  the  choice  of  some  random  probability  measure^ 
In  particular  we  consider  this  problem  when  C  is  a  disk  in  centered  at  0  and 

is  restricted  to  be  bivariate  normal.  Details  concerning  the  implementation  of  th 
estimation  procedures  and  a  Monte  Carlo  study  are  discussed  for  this  case.  This 
particular  problem  arises  when  we  are  concerned  with  assessing  the  accuracy  of  a 
targeting  procedure. 
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